A framework for solvation in quantum Monte Carlo 
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Employing a classical density-functional description of liquid environments, we introduce a rigor- 
ous method for the diffusion quantum Monte Carlo calculation of free energies and thermodynamic 
averages of solvated systems that requires neither thermodynamic sampling nor explicit solvent elec- 
trons. We find that this method yields promising results and small convergence errors for a set of test 
molecules. It is implemented readily and is applicable to a range of challenges in condensed matter, 
including the study of transition states of molecular and surface reactions in liquid environments. 



The physics of solvation, though poorly understood, 
plays a critical role in a wide range of systems from the 
biological to the technological. For example, the path- 
ways in protein folding [JJ-(3] and transition states of ionic 
reactions [H|S] are known to be highly solvent dependent. 
Development of a fundamental understanding of the ki- 
netics which underlie such processes requires an accurate 
description of the quantum mechanical processes involved 
in bond breaking and formation in solution. 

Unlike less rigorous electronic structure methods, 
quantum Monte Carlo methods [6] do provide the re- 
quired accuracy for transition states, reactants, and 
products needed to give reliable information on reaction 
pathways. However, a full quantum Monte Carlo treat- 
ment, in principle, requires solving Schrodinger's equa- 
tion for all of the involved solvent molecules, a difficulty 
radically compounded by the need to sample the phase- 
space of all thermodynamically relevant configurations 
of the solvent. This Communication introduces a frame- 
work for the treatment of solvation in diffusion Monte 
Carlo which completely eliminates the need for explicit 
solvent electrons and such phase-space sampling, while 
remaining completely ab initio and exact in principle 
in the same sense in which density-functional theory 
meets these criteria. 

Previous attempts to model the effects of solvation in 
quantum Monte Carlo fall into two categories: (i) sim- 
ulation of the environment through molecular dynam- 
ics [7J [5], and (ii) introduction of a polarizable contin- 
uum [9l [10]. The former approach, in principle, cap- 
tures molecular-scale effects such as solvation shells and 
non-local dielectric response. However, molecular dy- 
namics with empirical potentials [7J, though benefiting 
from a simplified description of the solvent, depends on 
a highly parameterized potential to describe the interac- 
tions between molecules. Additionally, such calculations 
require phase-space sampling, which makes the calcula- 
tions costly — especially if individual electronic structure 



calculations for the solute are needed for each fluid con- 
figuration. Approaches have been proposed to mitigate 
this expense within ab initio molecular dynamics calcu- 
lations [5], but this approach also is extremely expensive 
due to the need to include of all of the solvent electrons 
and nuclei. 

In contrast, the polarizable continuum model (PCM) 
is a computationally inexpensive tool for approximating 
solvation energies without phase-space sampling. Un- 
fortunately, this model lacks theoretical justification for 
the treatment of water as a continuum on molecular 
length scales, and so does not constitute a truly ab ini- 
tio method. Indeed, the pioneering work on solvation 
studies within diffusion Monte Carlo electronic structure 
methods [10] rested on the ad hoc introduction of a polar- 
izable medium and required a set of spheres with radii de- 
termined empirically. A further, practical disadvantage 
of the aforementioned PCM approach is that it requires 
potentially costly statistical evaluation of solvent poten- 
tials, and so it has yet to yield meaningful comparisons 
between predicted solvation energies and experiment. 

Joint density-functional theory |llj circumvents both 
the inaccuracies inherent in the polarizable continuum 
model and empirical molecular dynamics and the expense 
associated with ab initio molecular dynamics. This work 
presents an integrated approach to quantum Monte Carlo 
calculations within this new, rigorous statistical treat- 
ment of the solvent, and introduces a variational theo- 
rem which abrogates much of the computational com- 
plexity associated with achieving full self-consistency be- 
tween the solute and the solvent. While our approach 
in theory requires no adjustable parameters, we employ 
a simplified, approximate version of the theory in this 
first demonstration that does require a single empirically 
adjusted parameter. 

This Communication begins with a brief review of joint 
density-functional theory and then describes a quantum 
Monte Carlo formulation of that theory, compares predic- 
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tions to solvation data, and shows that self-consistency to 
within chemical accuracy is achievable with the computa- 
tional cost of a single quantum Monte Carlo calculation. 

Joint density -functional theory and quantum Monte 
Carlo. Joint density-functional theory [HI H2] states 
that, in principle, the exact quantum and thermodynam- 
ically averaged electronic and nuclear densities and the 
exact free energy of a solvated system can be obtained by 
minimization of a universal free-energy functional, with- 
out the need to sample explicitly the phase space of all 
possible configurations of solvent molecules. This mini- 
mization is carried out over all realizable average solute 
electron densities n(r) and average solvent site densities 
N a (r). (In the present work, N a (r) refers to the single- 
particle densities of the oxygen nuclei and protons com- 
prising the solvent, water.) Note that n(r) here refers 
only to the electron density of the solute because the 
electron density of the solvent need not be considered ex- 
plicitly. Naturally, the underlying indistinguishability of 
electrons implies that there is no unique decomposition of 
the total electron density into solute and solvent contri- 
butions. Consequently, the solution to the minimization 
problem is highly degenerate with many choices for n(r) 
leading to the same, exact free energy. This complication 
notwithstanding, the free energy and N a (r) obtained at 
the minimum are meaningful and represent the exact free 
energy and solvent site distribution of the combined sys- 
tem at equilibrium. (See Ref. |12) for a fuller discussion.) 
Finally, we note that the exact universal functional con- 
veniently decomposes into a sum of three terms, 

A[n, {N a }} = AmcW + $iq[{^4] + AA[n, {N a }}, (1) 

where AjjkM is the electronic Hohenberg-Kohn free- 
energy functional for the explicit system in isolation, $i q 
is the free-energy functional for the liquid when in isola- 
tion, and AA represents the coupling between the explicit 
system and the solvent. Note that, essentially being the 
definition of AA[n, {N a }}, Eq. is exact. 

Although Eq. @ is exact in principle, the forms of 
these three functionals on the right-hand side are un- 
known and need to be approximated in practice. Gener- 
ally, such approximations will break the aforementioned 
degeneracy and select a specific n(r) at the minimum of 
the approximated functional. In this Communication, we 
use diffusion Monte Carlo to describe the term AhkMi 
and, for our quantum Monte Carlo implementation of 
joint density-functional theory, we write the terms re- 
lated to the environment as 

A cnv [n] = mm (^[{N a }] + AA[n, {N a }}). (2) 

The variational derivative of Eq. ([I]) with respect to 
the exact electron density then yields the Euler-Lagrange 
equation = SA/Sn — S(Ahk + A eav )/dn, which is pre- 
cisely the equation for the isolated solute system in an 



external potential V cnv = 5A cnv [n]/5n. When an exter- 
nal potential is found for which the electron density n(r) 
apportioned to the solute yields back the same potential 
through this definition, the self-consistent, exact ther- 
modynamic state of the system will have been found. In 
principle, this exact solution can be obtained through 
multiple self-consistent iterations. The main difficulty of 
a quantum Monte Carlo implementation of Ahk in this 
process is the presence of statistical noise in the resulting 
electron densities, particularly in regions of low electron 
density, a problem which also plagues the approach of 
Ref. [TO]. There has been recent progress in reducing 
this noise [13] . but the results are not yet sufficiently 
clean for use in self-consistent calculations. We therefore 
now introduce a method to estimate the self-consistent 
solution of our solvation theory without evaluation of the 
electron density within quantum Monte Carlo. 

First-order estimator of self- consistency. A natural 
starting point for a joint density-functional theory quan- 
tum Monte Carlo calculation is a fully self-consistent 
joint density-functional theory calculation within the lo- 
cal density (LDA) or generalized gradient (GGA) approx- 
imation for Ahk- Such calculations provide trial wave- 
functions for quantum Monte Carlo calculations, and the 
resulting densities udft give estimates for the environ- 
ment potential V cnv — 8A env [n-£,YT\/5n. One can then 
perform a single quantum Monte Carlo calculation of the 
solute within this estimated potential, yielding both an 
initial quantum Monte Carlo density iiqmc and energy 
Aqmc[™qmc]- Note that the latter quantity is defined 
as the quantum Monte Carlo estimate of the functional 
Ahk', i-e., the energy of the solute without the energy 
associated with V env . 

The two calculations described in the above para- 
graph could then be combined to give a zeroth-order 
estimate of the free energy of the solvated system as 
Aqmc [ n QMc] + A cn v [^dft] ■ This, however, loses the ben- 
efits of the variational principle because the two terms are 
not self-consistent in that they are evaluated at different 
electron densities, and so the resulting estimate is not 
necessarily an upper bound for the final, converged re- 
sult. Also, the resulting error is first-order in the errors 
in the density, rather than second-order, as is usually as- 
sociated with variational calculations. We can correct for 
this by evaluating, with errors only in the second-order, 
A e nv[nQMc] = A env [tidft] + / d 3 rV env (uqmc - «-dft) 
because, by definition, V^ nv = SA eav [nDFT]/dn. The final 
result equals 

A = Aqmc [«qmc] + A cnv [?i D ft] + (3) 

J d 3 rV cnv (r) (iiQMc(r) - n DFT {r)) + 0(Sn 2 ), 

with errors that are second-order in both the difference 
between the exact converged density and the density- 
functional theory density «dft and the difference be- 
tween the exact converged density and the first-iteration 
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density ^qmc- Below, we show that this new approach 
is operationally nearly as good as full self-consistency, 
but with the effort of only a single quantum Monte Carlo 
calculation. 

Implementation. In this work we demonstrate our sol- 
vation technique with a simplified description of the en- 
vironment, but we stress that it may be employed just as 
readily with whatever functionals for the liquid become 
available. We note that the quality of our results now 
depend on our approximation for A env , which we hope 
to improve in the future. Specifically, we here employ the 
isodensity model of Petrosyan et al. [llj , which is similar 
to the dielectric model by Fattebert and Gygi [H] . These 
models take the dielectric constant to be a local func- 
tion of the electron density that switches smoothly from 
the value for vacuum at high solute electron densities to 
the dielectric constant of the bulk liquid for low elec- 
tron densities. This simplified joint density-functional is 
similar in spirit to the smoothly transitioning nonlinear 
polarizable continuum model described in Ref. [9], but 
here we work in a rigorous framework and also effectively 
achieve solute-solvent self-consistency. We have imple- 
mented our method in the open source code JDFTx |17j , 
interfaced with CASINO [T^]. In anticipation of ap- 
plying our method to large solvated surfaces, we use 
for our starting point the computationally inexpensive 
local-density approximation at full solute-solvent self- 
consistency. Fig. [ija) illustrates the cavity formation 
for acetone, a small molecule chosen for its well-known 
experimental solvation energy. Our reported solvation 
energies include, in addition to the electrostatic com- 
ponents from our simplified electrostatic functional, the 
cavitation energy shown in Table [I] estimated from clas- 
sical density functional theory following the procedure 
in Ref. [15] as included in the JDFTx package [17] and 
employing a functional based on Ref. [18] . 

The density-functional theory and Hartree-Fock cal- 
culations employ pseudopotentials from Burkatzki et al. 
[16) and expand the wave functions in a plane-wave ba- 
sis with a cutoff energy of 30 H (hartree). The local 
density approximation is employed for the density func- 
tional theory calculations, performed using JDFTx [17] . 
and a simulation box of 40 bohr 3 . The quantum Monte 
Carlo calculations are performed using CASINO jTHj and 
employ a trial wavefunction of a product of a single 
Slater determinant of density functional orbitals and a 
Jastrow correlation factor composed of electron-nucleus 
and electron-electron terms with expansion order eight, 
and electron-electron-nucleus terms with expansion order 
three, as described in Ref. [50] • The orbitals and exter- 
nal potential V^ nv are represented by B-splines. The pa- 
rameters of the Jastrow factor are optimized by variance 
minimization [21] . The diffusion Monte Carlo time step 
is 0.01 H _1 . (Going from 0.01 to 0.001 H _1 , we found the 
time-step error to be within the statistical uncertainty we 
report for the solvation energy for acetone in Figure [T] 
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FIG. 1. (Color Online) (a) Joint density-functional theory 
description of acetone in aqueous solution: electron den- 
sity contours as green (gray) surfaces, solvent as solid pur- 
ple (dark gray), (b) Solvation energies for acetone: diffu- 
sion Monte Carlo (DMC) as red (gray) x's, variational Monte 
Carlo (VMC) as green (gray) strikethrough x's, Joint Density 
Functional Theory (JDFT) as blue (gray) circles, Joint Den- 
sity Functional Theory best fit as a dashed black line, and 
experiment as a horizontal red (gray) line. 
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FIG. 2. (Color Online) Theoretical versus experimental sol- 
vation energies for ethanol, dimethyl ether (DME), propane 
(left to right): Joint Density Functional Theory as blue 
(gray) circles, zeroth-order diffusion Monte Carlo as black 
strikethrough x's, first-order corrected diffusion Monte Carlo 
as red (gray) x's. Experimental values from [23] . 



with the solvation parameter n c = 7 x 10~ 4 bohr - .) 
Molecular geometries are from the Computational Chem- 
istry Comparison and Benchmark Database: either from 
experimental data if available, or density-functional op- 
timization using the B3LYP functional and the cc-pVTZ 
basis set [2"2"] . 

For the continuum description of water, we use 
the local dielectric function e(n(r)) = 1 + (eb — 
l)erfc (ln(^)/( v / 2o-)^) /2, where e b is the bulk di- 
electric constant of the fluid, n c specifies the value of 
the solute electron density at which the dielectric cav- 
ity forms, and a controls the transition width [llj . 
This description of water leads to the expression 
A env [tt-dft] = §/ dr 3 (n DF T-^V)(-47r(V-e(n DFT )V)- 1 + 
47r(V 2 )~ 1 )(nDFT — N), where N represents the nuclei of 
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TABLE I. Formation energies (mH) for cavities at various 
values of the electron isodensity-contour parameter n c from 
the model of Ref. [11] , based on the cavitation energy method 
of Ref. [IS]. 



the solute. 

Results. Figure [ijb) shows the resulting solvation en- 
ergies for acetone as a function of n c . The results are 
not sensitive to er, so we leave this value fixed at 0.6 
and optimize the values of n c for use in diffusion Monte 
Carlo. All of the variational Monte Carlo results and 
nearly all the diffusion Monte Carlo results lie below the 
density-functional theory data. The variational Monte 
Carlo results with no Jastrow factor (and thus no cor- 
relation, but exact exchange energy) and the diffusion 
Monte Carlo results lie very near to each other, indicat- 
ing that the primary corrections in solvation energy to 
the density-functional results come from the exact treat- 
ment of the exchange and that corrections to correlation 
beyond the local-density approximation largely cancel for 
the solvation energies, at least for acetone. 

A least squares fit of the data shown in Fig. [TJd yields 
n c = 7.0 x 10~ 4 bohr~ 3 as optimal for diffusion Monte 
Carlo and n c = 8.1 x 10 -4 bohr -3 as optimal for use 
with density-functional theory. Figure [2] shows the re- 
sulting solvation energies of three molecules for both the 
zcroth-order expression and the first-order corrected ver- 
sion from Eq. |3]). The agreement between the quan- 
tum Monte Carlo results and experiment is encouraging 
given the particularly simple model employed here for 
the fluid. The figure also demonstrates the importance 
of using Eq. ^ to include the effects of self-consistency, 
particularly for ethanol. 

To estimate the remaining error between the corrected 
formula (Eq. ([3])) and full self-consistency, we employ two 
different electronic structure methods for which achieve- 
ment of full self-consistency is feasible and whose differ- 
ence in densities we expect to be similar to the density 
difference between density-functional theory and quan- 
tum Monte Carlo. We begin with an environment po- 
tential V env from a solvated density-functional theory 
(Hartree-Fock) calculation, and include it in a Hartree 
Fock (density-functional theory) calculation, attempt- 
ing to predict the final self-consistent energy using our 
proposed methods. Fig. [3] compares the zeroth and 
first-order corrected approximations with the fully self- 
consistent results when working this procedure in both 
directions. 

The data exhibit a number of behaviors which we ex- 
pect to be general trends. First, the first-order corrected 
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FIG. 3. (Color Online) Zeroth order approximation errors 
(squares) and first order approximation errors from Eq. Q 
(circles) relative to the exact self-consistent fluid-electronic 
structure minimization for Hartree-Fock calculations with 
LDA potentials in black (unfilled) and LDA calculations with 
Hartree-Fock potentials in red (gray) (filled), for the elec- 
trostatic solvation energies of ethanol, dimethyl ether, and 
propane (from left to right). 



data lie above the fully self-consistent result, regardless 
of the starting point. Second, the remaining (second- 
order and higher) errors are all quite small (0.2 mH or 
less) and on the order of about one-third the size of the 
first-order correction. Third, the remaining errors in the 
corrected results for a given molecule when going in either 
direction between density-functional theory and Hartree 
Fock are nearly identical. The difference between these 
remaining errors gives an estimate of the correction at 
odd orders (third and higher), indicating that the er- 
rors remaining are dominated by the second-order term. 
In aqueous solution, electrostatic screening, a negative- 
definite quantity, dominates this second-order term, thus 
explaining why the corrected results generally lie above 
the self-consistent solution. 

These observations strongly suggest that the self- 
consistency error in our corrected diffusion Monte Carlo 
results is less than one mH for ethanol, and even less for 
the other molecules. Our procedure thus likely gives an 
upper bound well within chemical accuracy of the results 
of full self-consistency, with the need for only a single 
quantum Monte Carlo calculation. 

Conclusion. The framework of joint density-functional 
theory provides a rigorous and efficient method for sol- 
vation in quantum Monte Carlo, without the need for 
phase-space sampling of the fluid. Moreover, a special 
procedure allows self-consistency of the joint calculation 
to be obtained (along with an estimate of the remain- 
ing errors) to within chemical accuracy, all at the cost 
of only a single quantum Monte Carlo total-energy cal- 
culation carried out in a fixed external potential. This 
solvation method applies just as readily to molecular and 
surface calculations. The procedure is general and may 
be used not only with quantum Monte Carlo but other 
correlated total-energy methods as well. 
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